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Abstract. Much of the work on flow through porous media, especially with regard to studies on 
the flow of oil, are based on "Darcy's law" or modifications to it such as Darcy-Forchheimer or 
Brinkman models. While many theoretical and numerical studies concerning flow through porous 
media have taken into account the inhomogeneity and anisotropy of the porous solid, they have 
not taken into account the fact that the viscosity of the fluid and drag coefficient could depend 
on the pressure in applications such as enhanced oil recovery. Experiments clearly indicate that 
the viscosity varies exponentially with respect to the pressure and the viscosity can change, in 
some applications, by several orders of magnitude. The fact that the viscosity depends on pressure 
immediately implies that the "drag coefficient" would also depend on the pressure. 

In this paper we consider modifications to Darcy's equation wherein the drag coefficient is a func- 
tion of pressure, which is a realistic model for technological applications like enhanced oil recovery 
and geological carbon sequestration. We first outline the approximations behind Darcy's equation 
and the modifications that we propose to Darcy's equation, and derive the governing equations 
through a systematic approach using mixture theory. We then propose a stabilized mixed finite 
element formulation for the modified Darcy's equation. To solve the resulting nonlinear equations 
we present a solution procedure based on the consistent Newton-Raphson method. We solve rep- 
resentative test problems to illustrate the performance of the proposed stabilized formulation. One 
of the objectives of this paper is also to show that the dependence of viscosity on the pressure can 
have a significant effect both on the qualitative and quantitative nature of the solution. 



1. INTRODUCTION 

Flow through porous media is commonly modeled using Darcy's equation [T]. Because of its 
popularity (in many fields like civil, geotechnical, petroleum engineering, composite manufacturing) 
the equation is incorrectly considered as a physical law, and is commonly referred to as "Darcy's 
law." However, it is important to note that Darcy's law is not a new balance law. It is an 
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approximation of the balance of linear momentum, and is valid under a plethora of assumptions. 
Some of the serious drawbacks of Darcy's law are 

• it cannot predict the stresses and strains in the solid, 

• it cannot predict the "swelling" and the attendant change in the pore structure, and 

• it cannot predict the induced inhomogeneity or anisotropy due to the deformation of the 



In fact, Darcy's law merely predicts the flux of the fluid through the porous solid. More importantly, 
this flux prediction is not accurate at high pressures and pressure gradients, which is one of the 
main themes of this paper. Henceforth, we shall use the terminology "Darcy's equation" instead 
of "Darcy's law." 

Over the years people have used Darcy's equation beyond its range of applicability. For example, 
Darcy's equation is used to model enhanced oil recovery and geological carbon sequestration. But 
in both these technologically important applications one has to deal with a high pressure range of 
10 MPa - 100 MPa. There is overwhelming evidence that viscosity is not constant and changes 
drastically with pressure in this range. 

Stokes [2] himself recognized that viscosity, in general, depends on pressure but that it can 
be assumed to be constant for a certain class of flows (such as pipe flows that involve moderate 
pressures and pressure gradients). Both enhanced oil recovery and geological carbon sequestration 
do not fall under this class of flows. Barus [3] suggested the following exponential relationship 
between viscosity fx and pressure p 



where /3 has units Pa . Andrade [I] proposed that the variation of viscosity with respect to density 
(p), pressure (jp) and temperature (8) can be modeled as 



where A, r and s are empirical constants. A monumental work that discusses in detail the effect of 
pressure on viscosity is the book by Bridgman [3] , which gives a comprehensive account of research 
concerning the physics of high pressures done prior to 1931. It is very important to realize that 
the flow characteristics of a fluid whose viscosity depends on pressure can be significantly different 



solid. 



(1) 



Kp) = W)exp[/3p] 



(2) 
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from that of the flow characteristics of a fluid with constant viscosity, and this forms one of the 
main subject matters of this paper. 

We shall now use Barus' formula to get a rough estimate of the variation of the viscosity (and 
hence the drag coefficient) with respect to pressure for some common organic liquids. For Naph- 
thenic mineral oil f) has been determined experimentally to be 23.4 GPa" 1 at 40°C [5]. Based on 
the Barus' formula for Naphthenic mineral oil we then have 

«t\ MP = 100 MPa)-Mp = 0.1MPa) 

(3) — — — — . x 100 = 935.7% 

w p(p = 0.1 MPa) 

It is worth remarking that since the viscosity varies exponentially, increasing the pressure range 
to that applicable to Elastohydrodynamics would lead to a value of the order of ten to the power 
of eight! On the other hand, using the Dowson-Higginson empirical formula [6], the density varies 
with respect to pressure as 

0.6p 



(4) p(p) = po 



1 + 



1 + 1.7p. 

where p is in GPa. Using the Dawson-Higginson formula we have 

p(p = 100 MPa) - p(p = 0.1 MPa) _ 
(5) p(p = 0.1MPa) " 5 - 1/0 

Many other experiments have also indicated that change in density due to changes in pressure at 
high pressures is indeed negligible. Thus, it would seem reasonable to neglect the effect of pressure 
on density and consider only the variation of drag function with respect to pressure. 

Despite experimental evidence of viscosity depending on pressure in an exponential manner, 
none of the studies of technologically important applications like enhanced oil recovery (EOR) and 
Carbon sequestration take such a dependence on pressure into account. Enhanced oil recovery 
(which is also referred to as tertiary or improved oil recovery) is employed after primary and 
secondary methods of oil recovery, which amount to only 20-40% of oil recovery. Using EOR 
techniques an additional 30% of oil can be obtained. One of the popular methods of enhanced oil 
recovery is gas injection. In the gas injection method, gas (typically, Carbon dioxide) is pumped 
into injection wells at high pressures, which pushes the oil to the surface at ejection wells. A 
typical enhanced oil recovery process is depicted in Figure [TJ In a subsequent section, using a 
typical reservoir simulation, we will show that the results (pressure, pressure gradients, and total 
flow rate) obtained by taking into account the exponential dependence of viscosity on pressure 
are qualitatively and quantitatively different from the classical results. To this end, we consider 
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modifying Darcy's equation to take into account the effect of pressure on viscosity. The resulting 
boundary value problem is solved using a stabilized mixed method. 

1.1. Stabilized mixed finite element formulations. Numerical methods for Darcy-type equa- 
tions can be classified into two main categories - primal (or single- field) formulations, and mixed (or 
two-field) formulations. In a primal formulation, the whole problem is written in terms of the pres- 
sure (which is considered as the primary variable). The governing equation will then be a Poisson's 
equation in terms of the pressure, and one can employ any of the standard numerical methods (for 
example, the Galerkin finite element formulation) to solve the resulting equation. Once the pres- 
sure field has been calculated, the velocity can be calculated using Darcy's equation by using the 
gradient of the pressure field. The two main disadvantages of using a primal formulation have been 
thoroughly discussed in the literature (for example, see [8]). Firstly, for low-order finite elements 
(which are often employed for generating complex meshes) the velocity is poorly approximated as 
one has to take the gradient of the pressure to calculate the velocity. In many situations, velocity is 
of primary interest. Secondly, primal formulations do not possess local mass conservation property 
with respect to the original computational grid. 

On the other hand, mixed formulations can alleviate the drawbacks of primal formulations, and 
are widely used in many numerical simulations. Some of the earlier works on mixed methods 
applied to flows through porous media are |2l EH IHl 021 Q31 OS 151 US1 QS1 H3 UB] - However, it is 
well-known that care should be taken when dealing with mixed formulations. A mixed formulation 
should meet the inp-sup (also known as LBB) stability condition \19\ 120] . This falls in the realm 
of stabilized formulations. A huge volume of literature is available on stabilized formulations for 
various equations (e.g., Darcy's equation, Stokes' flow, incompressible Navier-Stokes). A thorough 
discussion of stabilized methods is beyond the scope of this paper. An interested reader should 
refer to [2U 12^1 1251 IH1 HZl ISI] and. references therein, and to the following 
texts concerning stabilized methods [32j [33] . 

1.2. Main contributions of this paper. Some of the main contributions of this paper are 

• develop a new stabilized mixed finite element formulation for the modified Darcy's equation 
wherein the drag coefficient (which depends on the viscosity) is a function of pressure, 

• document that the dependence of viscosity on the pressure could have a significant effect 
on the nature of the solution (both velocity and pressure fields), and 
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• stress the need and importance of better models than Darcy's equation for modeling techno- 
logically important problems like enhanced oil recovery and Carbon sequestration for which 
Darcy's equation is physically not appropriate. 

1.3. An outline of the paper. The remainder of the paper is organized as follows. In Section[2]we 
derive Darcy's equation and our modification of Darcy's equation using mixture theory. In Section [3] 
we clearly outline our modifications to Darcy's equation. In Section [5] we present a stabilized mixed 
weak formulation for our modification of Darcy's equation, and also present a numerical solution 
procedure based on a consistent Newton-Raphson strategy. Representative numerical results are 
presented in Section [5j and conclusions are drawn in Section [6j 

2. ASSUMPTIONS BEHIND DARCY'S EQUATION AND MODIFIED DARCY'S 

EQUATION 

We now derive Darcy's equation using mixture theory (which is also referred to as the theory of 
interacting continua) [34 4 135 1 [36 ] I37|. This derivation reveals the main assumptions behind Darcy's 
approximation, and clearly demonstrates the range of its (in) applicability. Later, we will also derive 
a modification of Darcy's equation by relaxing one of the assumptions behind Darcy's equation. 
We first present the balance laws in mixture theory, and then consider the flow of a fluid through 
a (porous) solid. It should be noted that, unless explicitly stated, repeated indices do not imply 
summation. 

2.1. Derivation of modified Darcy equations. Consider a mixture of ./V constituents. Let 
and respectively, denote the density and velocity of the i-th constituent. The density of the 
mixture is defined as 

N 

(6) P = 5> (4) 

i=l 

The average velocity for the mixture (which is also referred to as the mixture velocity) is defined 
through 

(7) v = -J2p^v^ 

P i=l 

It should be noted that a variety of "mixture velocities" can be defined, and equation (|7|) is one 
of those possibilities. A detailed discussion about this important assumption and its implications 
with regard to the governing equations for mixture can be found in Rajagopal [38], and Rajagopal 
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and Tao [37]. Let denote the partial stress of the i-th constituent. The total stress is defined 
through 

N 



(8) ^ = E T 



t=l 



The total stress can also be defined in several ways, for a discussion of the same see Rajagopal and 
Tao [37] . The balance of mass can be written as 

(9) ^ + div[pW«W] = m« Vi 

where denotes the mass production of the i-th constituent. The balance of mass for the mixture 
as a whole warrants that 

N 

(10) ^ m « =0 

1=1 

The balance of linear momentum for the i-th constituent can be written as 

(11) (^^- + giad[v^]v^ = div[T« T ] + pW&W + I« 

where &® is the external body force acting on the i-th component, and I® denotes the interaction 
force that is exerted by other constituents on the i-th constituent in virtue of their being forced to 
co-occupy the domain of the mixture. By Newton's third law we have 

N 

(12) ^(lW+ m «„(9) =o 

i=l 

It is important to note that specification of J® is also part of the constitutive assumptions. In the 
absence of supply of angular momentum the balance of angular momentum can be written as 

N / N \ T 

(13) ^T«= HTT» 

i=i \i=\ ) 

which basically states that the total stress is symmetric. It should be noted that even in the absence 
of angular momentum supply the individual partial stresses need not be symmetric. 

We now consider the case of a mixture consisting of two constituents, namely, a fluid and a solid. 
The quantities associated with the solid and fluid have superscripts 's' and 'f, respectively. The 
assumptions behind Darcy's equation and their consequences are as follows: 
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(1) No mass production of individual constituents. That is, 

(14) m (s) = = 

(2) The solid is assumed to be a rigid porous media. Thus, one can ignore all balance laws for 
the solid. The stresses in the solid are what they need to be to ensure that the balance of 
linear momentum in the solid is met. 

(3) The flow is assumed to be steady 

(15) -^ = and — = 

(4) The fluid is assumed to be homogeneous and incompressible. Therefore, 

(16) grad[p (/) ] = and div[v^] = 

(5) The velocity and its gradient are assumed to be small so that the inertial effects can be 
neglected, which implies that 

(17) grad[t> (/) ]-u (/) = 

(6) The partial stress in the fluid is that for an Euler fluid (that is, viscous effects within the 
fluid are neglected), and the partial stress takes the form 

(18) TW = -p if) I 

(7) The only interaction force that comes into play is the frictional force at the boundaries 
of the pore. It is further assumed that this force is proportional to the relative velocity 
between the fluid and solid that is reflected by a drag-like term. Noting the fact that the 
solid is assumed to be rigid, and attaching the inertial frame to the solid, we have = 0. 
Hence, the interaction force on the fluid takes the form 

(19) = a (t>W - «W) = -avW 

The drag coefficient a, in general, can depend on pressure and relative velocity, v^ s > — v^'. 

(8) The drag coefficient a is independent of pressure and relative velocity, — 



Remark 1. In Darcy's equation the drag function is equal to the ratio of the viscosity of the fluid 
and the coefficient of permeability of the porous medium. That is, 

(20) a = ii 

where k is the coefficient of permeability. 

The above eight assumptions lead to Darcy's equation along with the incompressibility constraint 

(21a) av (/) + grad[p (/) ] = pb {f) 

(21b) div[t> (/) ] = 

Remark 2. One can obtain the Darcy-Forchheimer flow model |39| by relaxing the last assumption 
and allowing the drag coefficient to be a function of the magnitude of the relative velocity. Specifi- 
cally, a = ao + ai||t7^ ||, where a$ and oti are constants, and \\ ■ \\ is the 2-norm (or the Frobenius 
norm). Note that we have used the fact that = 0, and hence the magnitude of the relative 
velocity is equal to \\v^\\. 

In the next section, we shall relax one of the assumptions behind Darcy's equation by taking into 
account the dependence of viscosity on the pressure, and obtain a modification to Darcy's equation. 
We shall then show that even under this small (but realistic) extension, the response of the fluid is 
dramatically different to the the response of the fluid obtained using Darcy's equation. However, 
it should be notes that the framework offered by theory of interacting continua is quite general, 
and one can accommodate various aspects like mechanical deformation of the solid, thermal effects. 
In our future works, we shall show systematically how the response of various components in the 
mixture change by relaxing various other assumptions under Darcy's equation. 

3. GOVERNING EQUATIONS (MODIFIED DARCY EQUATION) 

The (spatial) position vector is denoted as x, and the gradient and divergence operators with 
respect to x are denoted as "grad" and "div", respectively. Let Vt C W ld be a bounded open 
domain, where "nd" denotes the number of spatial dimensions. The boundary of VL is denoted as 
r, which is assumed to be piecewise smooth. Mathematically, V is defined as V := Cl — $7, where 
fj is the closure of Q. We shall denote the velocity vector field as v{x). The pressure and density 
fields are denoted as p(x) and p(x), respectively. 
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As usual, r is divided into two parts, denoted by T D and T N , such that T D n T N = and 
r D U r N = r. r N is the part of the boundary on which normal component of the velocity is 
prescribed, and T D is part of the boundary on which pressure is prescribed. The modified Darcy 
equations can be written as 

(22a) a(p)v + grad[p] = pb in Q 

(22b) div[»] =0 in £1 

(22c) v{x) ■ n(x) = v n (x) on T N 

(22d) p(x)=po(x) onr D 

where a(p) (which has dimension of [ML _3 T -1 ]) is the drag function, po(x) is the prescribed 
pressure, v n (x) is the prescribed normal component of the velocity, b{x) is the specific body force, 
and n{x) is the unit outward normal vector to T. 

Remark 3. If T N = T then for well-posedness of (|22p one needs to satisfy the following compati- 
bility condition 



(23) / v n dr = 

Jr N =r 

which is a direct consequence of the divergence theorem. Also in this case, the solution for pressure 
is not unique. For uniqueness of the solution one needs to impose an additional constraint on the 
pressure. For example, 

(24) / p dn = 

Jn 

In numerical simulations, it is common to prescribe the pressure at a point, which is computationally 
easier than enforcing the constraint (|24l) , 



In this paper we consider the following three different drag functions 
(25a) a(p) = ao 

(25b) a(p) = a (1 + 0p) 

(25c) a(p) = ao exp | 



where ao > and j3 > are constants independent of p, v and p but can be functions of x. 
The dimension of the coefficient /3 is [M -1 LT -2 ]. The constant drag function (|25ap gives rise to 
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(standard) Darcy's equation, and the drag function (|25b|) is based on Barus' formula ([T]). Note 
that as p — y oo the drag coefficient ([25ap — > 0, drag coefficient (I25bj) — >• ao/3 (a constant), 
and drag coefficient (|25cp -4 oo. 

Several rigorous mathematical studies have addressed the equations governing the flows of fluids 
with pressure dependent viscosity. Existence of solutions have been established recently by Malek 
et al. |4Uj . Hron et al. [31], Franta et al. [32], and Bulicek et al. |43| . However, it is important 
to note that all the aforementioned existence results have been established under the condition 
^ —¥ constant as p — >■ oo. But, experiments suggest that ^ — > oo as p — > oo (for example, Barus' 
formula). Thus, rigorous results such as existence of solutions for fluids with pressure dependent 
viscosity are open with regard to the type of variation of the viscosity with pressure observed in 
experiments. 

It is, in general, not possible to obtain analytical solutions for the above boundary value problem 
(|22p . One has to resort to numerical solutions to solve realistic problems on complex domains. 
In this paper we employ the Finite Element Method (FEM), which is a powerful and systematic 
technique for numerically solving partial differential equations. To this end, we present the classical 
mixed formulation, and then present a new stabilized mixed formulation for solving modified Darcy 
equations. 

3.1. Classical mixed formulation. Let L 2 (d) be the space of square integrable (scalar) functions 
on ft, and h 2 (fl) be space of square integrable vector fields defined on O. For convenience, we denote 
the 1? inner-product over a spatial domain, K, as 

(26a) (a; b) K = [ a-b dK Va, b G L 2 (K) 

Jk 

(26b) (a;b) K = [ a ■ b dK Va,6eL 2 (iT) 

Jk 

The subscript K will be dropped if K is the whole of fl, that is K = CI. The (Hilbert) spaces 
H 1 ^) and H(div, O) are, respectively, defined as 

(27a) H\n) := {a(x) G L 2 {Cl) \ grad[a] G L 2 (f))} 

(27b) H(div,Q) := {a(x) G L 2 (fl) \ div[a] G L 2 (fl)} 
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Furthermore, we shall define the following function spaces 

(28a) V := {v(x) G H(div,0) | trace [v(x] ■ n(x)] = v n (x) on T N } 

(28b) W := {w(x) G H(div, 0) | trace [w(x) • n(x)] = on T N } 

(28c) V = L 2 (n), Q = il 1 (0) 

where trace[-] is the standard trace operator [20J. 

The classical mixed formulation for the modified Darcy equation (|22p can be written as: Find 
v(x) G V and p(x) G V such that 

(29) (w;a(p)v) — (div[w];p) + (w ■ n;p )r D ~~ (<7jdiv[V]) = (w; pb) \/w(x) G W, (/(a;) G P 

Assuming the domain and boundary are sufficiently smooth, for the case of constant drag function 
the above weak formulation (j29[) is well-posed |20j . A corresponding finite element formulation 
can be obtained by choosing suitable approximating finite element spaces, which (for a conforming 
formulation) will be finite dimensional subspaces of the underlying function spaces of the weak 
formulation. Let the finite element function spaces for the velocity, the weighting function associated 
with the velocity, and the pressure be denoted by V h C V, W h C W, and V h C V, respectively. A 
conforming finite element formulation of the classical mixed formulation reads: Find v h (x) G V h 
and p h (x) G V h such that 

(30) 

(w h ; a(p h )v h ) - (div[w h ];p h ) + (w h ■ n; p ) rD - (q h ; dw[v h }) = (w h ; pb) Vw h (x) G W h , q h (x) G V h 

For mixed formulations, the inclusions V h C V, W h C W, and V h C V are themselves not 
sufficient to produce stable results, and additional conditions must be met by the members of these 
finite element spaces to obtain meaningful numerical results. A systematic study of these types 
of conditions on function spaces to obtain stable numerical results is the main theme of mixed 
finite elements. One of the main conditions to be met is the Ladyzhenskaya-Babuska-Brezzi (LBB) 
inf-sup stability condition. For further details, see references |20|, 144] . 

It is well-known that under the classical mixed formulation many combinations of interpolations 
for velocity and pressure produce spurious oscillations (in the pressure field) when applied to prob- 
lems involving incompressibility as a constraint \20\ 144] . In particular, the equal-order interpolation 

for the velocity and pressure (which is computationally the most convenient) does not satisfy the 
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LBB condition, and produces spurious oscillations in the pressure field. However, elegant solutions 
for approximating the velocity and pressure that are stable under the classical mixed formulation 
have been proposed [1SJ [Ml H7J H5J [Ml [SD]. These discrete spaces have been successfully used in 
many numerical simulations, and good accuracy has been attained in both the velocity and pres- 
sure fields. But a computer implementation of these methods is more involved because of different 
interpolations for the velocity and pressure. 

In the next section we present a new stabilized mixed formulation for modified Darcy equations 
under which the computationally convenient equal-order interpolation for velocity and pressure is 
stable. 



A variational multiscale mixed formulation has been proposed and studied by Hughes and Masud 
[8] and Nakshatrala et al. [T7] for Darcy's equation. In reference [T7] it has been shown that this 
variational multiscale formulation for Darcy's equation passes three-dimensional patch tests even 
for distorted meshes, and performs well even on complex geometries with unstructured meshes. In 
this paper we extend this formulation to the modifications to Darcy's equation developed in this 
paper. 

A stabilized mixed formulation for modified Darcy's equation (|22p can be written as: Find 
v(x) G V and p(x) € Q such that 

(w;a(p)v) - (div[w];p) + (w ■ n;p ) rD - (q;dw[v}) - (w; pb) 



The terms in the first line of equation (|3ip are from the classical mixed formulation, and the terms 
in the second line are the stabilization terms. The factor 1/2 on the second line is the stabilization 
parameter, which is determined neither by mesh-dependent parameters nor on material properties 
(like viscosity). The stabilization term is based on the residual of the modified Darcy's equation 
(|22p . Similarly, one may add a residual-based term using the incompressibility constraint, which in 
some cases improves stability (e.g., see [8]). We now solve the above weak formulation pip using 
the Finite Element Method. 



4. A STABILIZED MIXED FORMULATION 



(31) 
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4.1. A finite element approximation. Let the domain be decomposed into "Nele" non- 
overlapping open element subdomains. That is, 

Nele 

(32) n = \J n e 

e=l 

where a superposed bar denotes closure. The boundary of il e is denoted as <9il e := Cl e — £l e . For 
a non-negative integer m, T >m {Vt e ) denotes the linear vector space spanned by polynomials up to 
mth order defined on the subdomain Q e . We shall define the following finite dimensional vector 
spaces of V, W and Q: 

(33a) V h := iv h (x)eV\ v h (x) G (C°(fi)) nd , v h {x) | Qe G (V*^ 6 ))^ , e = 1, ■ ■ ■ ,WeZej 

(33b) W h := |w ft (aj) G W | io h (aj) G (C°(0)) nd , w h {x)\^ G (^(O 6 ))™* , e = 1, • • • ,iVeZej 

(33c) Q h := [p h (x) G Q \ p h (x) G C ^),/^)^ G V l (n e ), e = 1, • • ■ ,iVe/e} 

where A; and / are non-negative integers, and recall that "nd" denotes the number of spatial di- 
mensions. A corresponding finite element formulation can be written as: Find v h {x) G V h and 
p h (x) G Q h such that 

(w h ;a(p h )v h ) - (dw[w h ];p h ) + (w h -n;p Q ) r u - (q h ; div[v h ]) - (w h ;pb) 
(34) -i (a(/)w /l + grad[g h ];Q- 1 (/)(a(/)- i ; /l + grad[/] - p6)) =0 Vw h G W' 1 , g h G Q h 

Due to the presence of the term |(grad[g^]; a~ 1 (p h )grad[p h ]), the above formulation circumvents 
the LBB condition, and hence one can employ equal-order interpolation for velocity and pressure, 
which will be confirmed by numerical results presented in a subsequent section. 

We solve the resulting nonlinear equations based on a consistent Newton- Raphson solution pro- 
cedure. We now derive corresponding residual vector and tangent matrix, which are required for 
such a procedure. 

4.2. Residual vector. We shall define the nodal unknowns for the velocity (denoted as v e ) and 
pressure (denoted as p e ) for a given element f2 e as 







v\,\ ■ 


■ v 1>nd 




Pi 


(35) 


Ve ■ = 






, Pe ■ = 








_ Vn,l ■ 


Vn,nd 
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where n is the number of nodes in the given element, and nd (as mentioned earlier) is the number 
of spatial dimensions. The solution, v(x) and p(x), over the element O e are interpolated in terms 
of corresponding nodal values (which have a superposed hat) as 



(36) 



v(x) = v^N T (x), p(x) = N(x)p e 



where N(x) is a row vector of shape functions. Likewise, the weighting functions, w(x) and q(x) 
over an element are interpolated as 



(37) 



w(x) = w^N T (x), q(x) = N(x)q e 



where the nodal (arbitrary) weights w e and q e are defined similar to v e and p e (which are defined in 
equation ()35|) ) . One can construct residual vectors corresponding to equation (|3ip by substituting 
equations (|36p into equation (|3ip . and invoking the arbitrariness of w e and q e . The element residual 
vectors can be written as 

R e v (v e ,p e ):= [ (N T I)a(p)v dfl — [ vec [B T ] p dfi + / (N T Q I)np dT 

(38) - / (N T I)pb dn~l [ (N T /) (a(p)v + grad[p] - pb) dQ 

iff 2 J U e 

(39) R e Jv e ,Pe) ■=- f iV T divM dn - I [ B a~ l {p) (a{p)v + gmd[p] - pb) dtl 

where v and p are approximated using equation (|36p . vec[-] is an operation that represents a matrix 
as a vector, is Kronecker product [51] (see Appendix [7J), B := DNJ~ , J := dx/d£ is the 
element Jacobian matrix, DN represents a matrix of (first) derivatives of element shape functions 
with respect to reference coordinates For example, for a three-node triangular element the matrix 
DN will read 



(40) 



DN 



dNi dNi 



dN 3 dNa 
L 96 dS,2 



We shall define the global unknown vectors v and p as 







vi,i ■ 


■ v 1>nd 




Pi 


(41) 


v := 






, p : = 










Vnn,nd 




Pnn 
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where nn denotes the total number of nodes in the finite element mesh. Global residual vectors 
can be constructed from element residual vectors using the standard assembly procedure 

ANele « Nele 

e=l K(Ve,P e ), Rp{v,P) = A e=1 R e p {v e ,p e ) 

where A is the standard assembly operator (Recall that Nele denotes the number of ele- 

ments.) The finite element solution can be obtained by solving simultaneously 



(43) 



Ry(v,p) = and R p (v,p) = 



4.3. Tangent matrix. Using a Newton-Raphson type approach, one can obtain the solution in 
an iterative fashion using the following update equation until the residual is under a prescribed 
tolerance 



(44a) 


vec 




= vec 




+ vec 




(44b) 


vec 




= vec 




+ vec 





The updates at iteration i are calculated from the following system of equations 
(45) 







1 


vec 






hi 


Rv } 




-Kpp 




vec 


Ap« 






{ Rp \ 



where 



ANele j» Nele 

e=l K lvi K vp = " e =l K%p, Kpv 

and the element matrices are defined as 



Nele 

'e=l "-pvi "-pp 



Nele 



A K 

* *e=l pp 



(47a) K% v = \ I (N T I) a(p) (N /) dft 

(47b) K e = i / (AT T 0j)« ( ^ ) Ndn~ [ vec [S T ] AT dft - 1 / (N T Q I) B T dtt 

(47c) K e = - [ AT T vec [S T ] T dfi - J / S (AT I) df) 

(47d) = -1 ^ Ba~\p)B T dSl + ±J B (grad[p] - p6) AT dfi 

Remark 4. From equation (|47p z'i is evident that for the modified Darcy's equation, in general, we 
have Kp V ^ K% V T , and the matrix K e pp is not symmetric. But in the case of Darcy's equation (for 
which a is independent of pressure) we do have K^ v = K^ p T , and the matrix K e vv is symmetric. 
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It should be noted that symmetry of tangent matrix (or the lack of it) is an important factor to be 
considered in the selection of solvers for the resulting system of linear equations, and also in the 
design/selection of preconditioners for iterative solvers, which are commonly employed in large-scale 
problems. 

5. REPRESENTATIVE NUMERICAL RESULTS 

In this section, we present representative numerical results to illustrate the performance of the 
proposed stabilized mixed formulation. In all our numerical simulations we have employed equal- 
order interpolation for pressure and velocity as shown in Figure [2j 

First, we shall non-dimensionalize the governing equations. To this end, we define the following 
non-dimensional quantities (which have a superposed bar) 



,. Q s x _ v _ p _ a _ a _ p - b - _ v n _ p 

(48) x = -, v = — , p = — , a = , a = , p = , b = — , j3 = @P, v n = — , p = 

L V P a vei a rei p r ef B VP 



where L, V, P, ct re f, p re f and B respectively denote reference length, velocity, pressure, drag 
coefficient, density and specific body force. The gradient and divergence operators with respect to 



x are denoted as "grad" and "div" , respectively. The scaled domain O sca i e d is defined as follows: a 
point in space with position vector x € ^ sca ied corresponds to the same point with position vector 
given by x = xL £ 0. Similarly, one can define the scaled boundaries: dr2 sca ied> r^. aled , and r^, aled . 
The above non-dimensionalization gives rise to two dimensionless parameters 



(49) A:=^^, C:=^ 



A corresponding non-dimensionalized form of the drag functions given in equation (j25[) can be 
written as 

(50a) a(p) = Qo 

(50b) a(p) = Q (1 + fa) 

(50c) a(p) = aoexp[/3p] 
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We shall write a non-dimensional form of the modified Darcy equation as 



(51a) A a(p)v + gradfp] =C pb in f2 sca i cd 

(51b) div[v) =0 in S7 sca i cd 

(51c) v-n = v n (x) onrf calcd 

(51d) p(x)=po(x) onrf calcd 

where fi(x) is the unit outward normal to the boundary <9f2 sca i e d- 

5.1. A simple one-dimensional problem. We shall test the proposed stabilized mixed formula- 
tion using a simple one-dimensional problem, which is pictorially described in Figure [3j Pressures 
of p\ and p2 are respectively prescribed at the left and right ends of the unit domain, and we neglect 
the body force. The governing equations for this test problem can be written as 

(52a) Aa(p)v(x) + ^ = in (0, 1) 

ax 

(52b) $^ = in (0,1) 

ax 

(52c) p(x = 0) = pi, p(x = 1) = p2 

For this test problem, the analytical solutions with the drag functions defined in equation (|25p can 
be written as 



(53a) for the case a(p) = ao < 



p(x) = (p 2 -pi)x+pi 

(P2-Pl) 



v \ x ) ~- MT 



(53b) 



for the case a{p) = ao(l + j3p) 



p(x) 



v[x 



{l + pplf *{l + pp 2 ) x -l 



Aa (3 



hi 



1+/3P2 
l+Ppl 



(53c) for the case a(p) = ao exp 



P( x ) = -J- m {(1 - x ) ex P [-fiPi] + a^exp [-/3p 2 ] } 
, = aE^ ( ex P [~?Pi[ ~ ex P [-PPi] } 

In Figure 2] we compare the numerical solutions against the analytical solutions for various drag 



functions. In all the considered cases, the proposed numerical formulation performed well. 
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5.2. Five-spot problem. This subsection presents numerical results for a quarter of the five- 
spot problem, which exhibits elliptic singularities near (injection and production) wells and is 
widely used as a good model problem to test the robustness of a numerical formulation. Taking 
into account the quarter symmetry in the five-spot problem, the source and sink strengths at 
injection and production wells are, respectively, taken as +1/4 and —1/4. There is no volumetric 
source/sink (i.e., b(x) = 0). The five-spot problem is pictorially described in Figure [SJ The 
numerically obtained pressure profiles using Darcy's equation and modified Darcy's equation with 
the viscosity given by Barus' formula are shown in the Figures [6] and [7J Firstly, there are no 
spurious oscillations in the pressure field even for the modified Darcy's equation, and the proposed 
mixed stabilized formulation performed well. Secondly, for this test problem, pressure in the case 
of the modified Darcy equation exhibits steeper gradients compared to the pressure field based on 
Darcy's equation. In Table Q] we have shown that the norm of the residual vector exhibits terminal 
quadratic convergence, which basically shows that the tangent matrix corresponding to the given 
residual vector is correctly formulated and calculated. In nonlinear problems, (due to the lack of 
analytical solutions) the terminal quadratic convergence is considered as a good measure to check 
the computer implementation of a numerical formulation. In Figure [8] we compared the number 
of Newton-Raphson iterations for various values of /3 using four-node quadrilateral elements for 
a relative tolerance of 10 -12 . In the same figure we have also shown the variation of exp[/3p max ] 
with respect to /3, where p max is the maximum pressure in the computational domain. (Note that 
Pmax occurs at the injection well.) As one can see, as j3 increases exp[/3p max ] increases drastically, 
and very steep pressure gradients occur near the injection and production wells. To obtain results 
for much higher values of /3 one may have to precondition the resulting linear equations at every 
Newton-Raphson iteration in order to avoid ill-conditioned matrices, and also employ a finer mesh. 

Remark 5. It is well-known that the Newton-Raphson solution procedure will not always exhibit 
terminal quadratic convergence. For example, if the multiplicity of the desired root is greater than 
one, the Newton-Raphson solution procedure may not exhibit terminal quadratic convergence. Also, 
in some cases, the Newton-Raphson procedure may not even converge. For a detailed discussion of 
the convergence properties of the Newton-Raphson solution procedure see Reference |53j . 

5.3. /i-convergence analysis for a two-dimensional problem. Let us consider a bi-unit square 
[0, 1] x [0, 1] as our computational domain. Let us assume that the velocity and pressure are 
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Table 1. Five-spot problem: Convergence of Newton-Raphson solution procedure 
for solving modified Darcy's equation using four-node quadrilateral (Q4) and three- 
node triangular (T3) elements. We have used Barus' formula with cxq = 1 and 
/3 = 0.3. The test problem is pictorially described in Figure [5j As one can see from 
the table, the solution procedure exhibits terminal quadratic convergence. 



Iteration # 


Norm of residual (Q4) 


Norm of residual (T3) 





0.242406260 


0.291915904 


1 


0.082035285 


0.086348093 


2 


0.037068940 


0.037721275 


3 


0.004513442 


0.004318753 


4 


4.1044706e-05 


3.4343232e-05 


5 


1.4244950e-09 


7.6123308e-10 


6 


2.0953922e-15 


3.4935047e-15 



respectively given as 
(54) 

v x (x, y) = sin(7rx) cos(iry), v y (x, y) = — cos(7rx) sm(ny), p(x, y) = 1 + 25xy(x — \)(y — 1) 

We shall assume that p = 1, c*o = 1, and /3 = 2. Using the above assumed solution, the required 
body force under the Barus' formula ([I]) should be equal to 

(55a) b x (x, y) = exp[/3(l + 25xy(x — l)(y — 1))] sin(7rx) cos(7r?/) + 25(2x — l)y(y — 1) 

(55b) b y (x, y) = — exp[/3(l + 25xy(x — l)(y — 1))] cos(-7ra;) sin(7ry) + 25x(x — l)(2y — 1) 

and the boundary conditions are given by 

(56) v x (x = 0,y)=v x (x = l,y) = Q, v y (x,y = 0) = v y (x,y = 1) = 

Using the aforementioned boundary value problem, we shall perform numerical convergence studies 
of the proposed stabilized mixed formulation for both four-node quadrilateral and three-node trian- 
gular elements. Typical uniform finite element meshes used in the /i-convergence analysis are shown 
in Figure [9j Note that all the elements in the four- node quadrilateral mesh are squares, and all the 
elements in the three-node triangular mesh are isosceles right-angled triangles. In the numerical 
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Table 2. Unstructured three-node triangular mesh: Convergence of Newton- 
Raphson solution procedure for solving modified Darcy's equation. We have used 
Barus' formula with ao = 1 and (3 = 2. As one can see from the table, the solution 
procedure exhibits terminal quadratic convergence. 



Iteration ft 


Norm of residual 





1.353515824 


1 


3.594631185 


2 


1.781220654 


3 


31.25838881 


4 


0.169839581 


5 


0.015076303 


6 


3.125925098e-05 


7 


7.498072058e-10 


8 


6.532312307e-16 



convergence studies presented in this subsection, we have taken h to be equal to the length of the 
side for square elements, and length of the base (or height) for the isosceles right-angled triangles. 

The contours of pressure and velocity are shown in Figures [TO] and [TU The rates of /i-convergence 
of these finite elements for the aforementioned problem are shown in Figure [12] (We have used the 
natural logarithm in these figures.) For the chosen problem, the proposed formulation converges 
in both pressure and velocity with respect to /i-refinement. The above problem (|55j) is also solved 
using an unstructured triangular mesh. The obtained numerical results and mesh are shown in 
Figure \13\ quadratic convergence of the Newton- Raphson solution procedure for this problem is 
illustrated in Table [2j and the proposed mixed stabilized formulation performed well. 

5.4. A typical reservoir simulation. Figure Q] shows how enhanced oil recovery is achieved in 
the field. A corresponding computational idealization of enhanced oil recovery is shown in Figure 
[T4"l which will be discretized using finite elements. The computational mesh is shown in Figure [T5l 
It should be noted that the corner points at the production well are reentrant corners, and hence 
the velocity v (which is proportional to gradient of pressure) is infinity at those reentrant corner 
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points. We now compare the solutions obtained using Darcy's equation (constant drag function) 
and modified Darcy's equation using Barus' formula. 

For this problem, we have taken So = 1, P = 0.005, P = 1 atm, B = 9.81 m/s 2 , b = [0, — 1] T , 
A = 1 and C = 1. The top surface at the production well is at atmospheric pressure (see Figure 
PT1| . which is given by the non-dimensional parameter p a t m — 1 ( as we have taken the reference 
pressure to be P = 1 atm). 

In Figures [16] and [T71 pressure contours are shown for the cases p e nh = 5 and p en h = 500, 
respectively. For low prescribed pressures at the injection wells, the pressure profiles obtained 
using Darcy's equation and modified Darcy's equation are similar. However, this is not the case for 
higher pressures at the injection wells, and the pressure profiles using Darcy's and modified Darcy's 
equations are significantly different. The pressure using modified Darcy's equation exhibit steeper 
gradients near the injection wells. It should be noted that pressure gradients play an important 
role in the study of fracture of porous medium. 

In Figure [T8l the variations of pressure with respect to the horizontal distance from the produc- 
tion well at various depths are plotted for Darcy's and modified Darcy's equations. In this study 
we have taken p en h = 1000. As one can see from Figure [TBI the pressures from Darcy's equation 
and modified Darcy's equation are qualitatively and quantitatively different. Most importantly, 
the cone shapes are completely different. (Cone shape is the variation of pressure with respect 
to the horizontal distance from the production well, and is considered an important parameter in 
assessing the well performance and also regulating reservoir-operating conditions.) For the case of 
Darcy's equation, the graph is concave upwards, while the graph using modified Darcy's equation 
(with Barus' formula) is convex upwards. 

In Figure [19] we compared the total flux at the production well for various pressures at the 
injection wells (which is given by p e nh) using Darcy's equation and modified Darcy's equation 
(based on Barus' formula). Even for this case the pressure at the production well is p a tm = 1- The 
total flux at the production well is calculated using 



where the part of boundary is indicated in Figure O As one can see in Figure \19\ modified 
Darcy's equation predicts a ceiling for the total flux with respect to pressure, which is what is 



(57) 
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expected physically. On the other hand, Darcy's equation predicts continued linear increase of the 
total flux with respect to the pressure. 

One can expect even more interesting departures from the classical results when the Darcy- 
Forchheimer and Darcy-Forchheimer-Brinkman equations are modified to incorporate the effect of 
pressure on the viscosity and drag function. Studying the interactions between the solid deformation 
and fluid flow, with the effects due to the pressure being taken into account will provide even more 
insights into the real problem, and is part of our future work. 



6. CONCLUSIONS 

In this paper we have studied the flow of an incompressible fluid through a rigid porous solid under 
high pressure and pressure gradients. For such problems, Darcy's equation is not a good model as it 
assumes constant viscosity and drag coefficient, which is contrary to experimental evidence. Here, 
we have considered a modification to Darcy's equation wherein the drag coefficient is a function of 
pressure, which is the case with many liquids for large range of pressures. We have allowed the drag 
coefficient to vary with pressure in different ways, which are primarily motivated by experimental 
observations. We have also presented a new stabilized mixed formulation for the modified Darcy's 
equation along with the incompressibility constraint. The proposed formulation allows equal-order 
interpolation for pressure and velocity, which is not stable under the classical mixed formulation. 
A noteworthy feature of the formulation is that the stabilization parameter neither involves mesh- 
dependent nor material parameters. We have presented representative numerical results to show 
that the proposed formulation performs well. We have also shown that the solution (both velocity 
and pressure fields) based on the modified Darcy's equation is significantly different from the 
solution based on Darcy's equation. In particular, using a representative reservoir simulation, we 
have shown that the Darcy's model may predict significantly higher flow rate at high pressure as 
it does not account for variation of viscosity with respect to pressure. 



7. APPENDIX 

7.1. Notation and definitions. We now define some quantities that will be useful in writing the 
finite element residual vector and tangent stiffness matrix in a compact manner. Let A and B be 
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matrices of size n x m and p x q, respectively. That is, 



ai, 



&1,1 • • • &l,q 



3p,l 



The Kronecker product of these matrices is an np x mg matrix, and is defined as 



AQB 



ai^-B . . . ai >m B 



a ni iB . . . a n ^ m B 

Note that the Kronecker product can defined for any two given matrices (irrespective of their 
dimensions). The vec[-] operator is defined as 



vec[A] := 



ai, r 



Some relevant properties of Kronecker product and vec[-] operator are as follows: 

• vec[ACB] = (B T A) vec[C] 

• (A B) (C D) = {AC BD) 

• vec[A + B] = vec[A] + vec[S] 

For a detailed discussion on Kronecker products and vec[-] operator see References [51 |, 154 1 [55] . 
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Figure 1. A schematic diagram of enhanced oil recovery. This picture is taken 
from https:/ /www. I Inl.gov/str/NovemberOl/Kirkenda 1 1, html. 




Figure 2. Typical two-dimensional finite elements employed in the numerical sim- 
ulations, and the location of velocity and pressure unknowns. Equal-order interpo- 
lation is used for the velocity and pressure. 

p(0) = pi p(l) = p 2 

>* 

1.0 



Figure 3. Schematic description of the one-dimensional problem. 
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a(p) = 1.0 

a(p) = 1.0 + 0.01 p 

a(p) = 1.0 + 0.05 p 

a(p) = exp(0.01 p) 

a(p) = exp(0.05 p) 

Analytical 



'Sh 100 




Figure 4. One-dimensional problem: The numerically obtained pressure profiles 
are compared against analytical solutions for various drag functions. We have chosen 
pi = 200, p2 = 1, .4. = 1, and ao = 1. The body force is neglected. The test problem 
is pictorially described in Figure El and the analytical solutions are given in equation 
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Figure 5. Pictorial description of a quarter of five-spot problem. 
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(a) Pressure contours using Darcy's equation 




(b) Pressure contours using modified Darcy's equation 

Figure 6. Five-spot problem using 400 four-node quadrilateral elements: This fig- 
ure compares elevation plots for the pressure profile using Darcy's equation (top) 
and modified Darcy's equation (bottom). The pressure based on modified Darcy's 
equation with the viscosity given by Bar us' formula exhibits steeper gradients com- 
pared to ones obtained using Darcy's equation. As one can see there are no spurious 
oscillations in the pressure field despite the steep gradients in pressure near the in- 
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jection and pressure wells. For Barus' formula we have used qq = 1 and /3 = 0.3. 




(a) Pressure contours using Darcy's equation 




(b) Pressure contous using modified Darcy's equation 

Figure 7. Five-spot problem using 800 three-node triangular elements: This figure 
compares elevation plots for the pressure profile using Darcy's equation and modified 
Darcy's equation (bottom). The pressure based on modified Darcy's equation with 
the viscosity given by Barus' formula exhibits steeper gradients compared to ones 
obtained using Darcy's equation. As one can see there are no spurious oscillations 
in the pressure field despite the steep gradients in pressure near the injection and 
production wells. For Barus' formula, we 3 have used ckq = 1 and j3 = 0.3. 
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Figure 8. Five-spot problem using 400 four-node quadrilateral elements: The top 
figure shows the number of Newton-Raphson iterations for various values of f3. 
The bottom figure shows the variation of exp[/5p max ] for various values of /3, where 
exp[/3p max ] denotes the maximum pressure in the computational domain. We have 
employed Barus' formula with ckq = 1. 
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(a) (b) 

Figure 9. Typical meshes used in the /i-convergence analysis: (a) four-node quadri- 
lateral (b) three- node triangular (right) meshes. 
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(a) Using four-node quadrilateral elements (b) Using three-node triangular elements 

Figure 10. Structured meshes: contours of pressure using 100 uniform four-node 
quadrilateral (left) and 200 three-node triangular (right) elements. There are no 
spurious oscillations in the pressure field. 
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(a) z-velocity using Q4 elements (b) a;-velocity using T3 elements 




(c) j/-velocity using Q4 elements (d) i/-velocity using T3 elements 



Figure 11. Structured meshes: contours of x- velocity (top) and y-velocity (bot- 
tom) using 100 uniform four-node quadrilateral (Q4) elements and 200 three-node 
triangular (T3) elements. 
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(a) Convergence for four-node quadrilateral elements 



—eL 

2.5 
-3 
3.5 
-4 
















-•-L 2 error of p 


















error of v 












lope = 1 


7 




























4.5 
-5 








s 


ope = 1 . 


7 

















































cl I I I I i I I I I I 

1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 



-log(h) 

(b) Convergence for three-node triangular elements 

Figure 12. 2D /i-convergence analysis: Rate of convergence of pressure and veloc- 
ity in the L 2 (tt) norm for four-node quadrilateral and three-node triangular finite 
elements. 
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(a) Computational mesh 
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(b) Pressure contours 




(c) x-velocity contours 



(d) y-velocity contours 



Figure 13. Unstructured three-node triangular mesh: mesh (top-left), pressure 
(top-right), x-velocity (bottom- left) and y-velocity (bottom-right). Barus' formula 
is employed with qq = 1 and /3 = 2. 



.35 



Injection well 

^ A 



Production well 



W 



Injection well 




it 



^ • n = o rp 
_^ / 



v ■ n = 



r N P (x) = p atm r 



N 



r D 

1 2 



~7 7" 



W = 0.1 

— 7 

L = 2 



r D 

1 2 



i> ■ n = 



II 

o 
B 



Figure 14. Top figure is a pictorial description of a reservoir. Water, steam or 
Carbon dioxide is pumped through injection wells, and crude oil is collected at the 
production well. Bottom figure shows the computational idealization of the reservoir 
with appropriate boundary conditions. On we prescribe p(x) = p a tm ; and on 
we prescribe p(x) = p e nh- 
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Figure 15. A four- node quadrilateral (structured) finite element mesh used in the 
reservoir simulation. There are 8000 elements, 8241 nodes, and each node has three 
degrees-of-freedom (x- and y- velocities, and pressure). 
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Figure 16. Pressure contours using Darcy's equation (top) and modified Darcy's 
equation (bottom) for p en ^ = 5. 
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Figure 18. Variation of pressure with respect to horizontal distance from the pro- 
duction well at various depths of the reservoir using Darcy's equation and modified 
Darcy's equation for p cn h = 1000. 
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FIGURE 19. Total flow rate at the production well using Darcy's equation and 
modified Darcy's equation (using Barus' formula) for various values of p e nh- 
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